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Abstract 

System matrices for Euler-Bernoulli beam problems for the meshless local 
Petrov-Galerkin (MLPG) method deteriorate as the number of nodes in the beam models 
are consistently increased. Tire reason for this behavior is explained. To overcome this 
difficulty and improve the accuracy of the solutions, a local coordinate approach for the 
evaluation of the generalized moving least squares shape functions and their derivatives 
is proposed. The proposed approach retains the accuracy of the MLPG methods. 
Introduction 

Meshless methods are increasingly being viewed as an alternative to the finite 
element method [1-3]. Recently, a meshless local Petrov-Galerkin (MLPG) method has 
been presented for C° and C 1 problems [3,4], hr these methods, moving least squares 
(MLS) interpolants [1] are used for C° problems and generalized MLS interpolants are 
used for C 1 problems [4], References 3 and 4 showed excellent performance of the 
MLPG method for potential and elasticity problems and a good performance for beam 
problems. 
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When all the chosen parameters in the MLPG method are held constant and the 
number of nodes in the models are consistently increased, the error norms do not 
decrease; rather they show increases compared to coarser idealizations. The reasons for 
this behavior are studied. A local coordinate approach to the MLS interpolation is 
proposed. The proposed local coordinate approach is implemented and evaluated by 
applying it to three simple test cases. 

Behavior of the MLPG method with mesh refinement 


Tire notation of reference 4 is used in this note for brevity and convenience in 
presentation. The MLPG equations are 

K (n° de) d + K^d - f/ n0de) - f, (My) =0 (1) 

where 

{d }={u, e\ u 2 , e 2 , ...} (2) 

are the fictitious nodal values of deflections u and slopes 6, and the matrices in Eq. (1) 
are defined as in Eq. 35-36(g) of reference 4. 

Tire MLPG equations are derived using a weighted residual weak form of the 
governing equations. Tire trial functions used for the beam problems are derived using 


the generalized MLS interpolation [4] as 
u{x) = ^u t yf\ u \x) + G t yr ( f\x) 

i = 1 

where 

¥i U) (x)= X py(*)[A- 1 P T wl /7 

(6) m f I t 

¥i ) (x)= l pj(x )[ A l Pj w{ /7 

with 

[A] = P T wP + PjwP x . 


( 3 ) 

( 4 ) 

( 5 ) 
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Ill Eq. (5) P is an (n,m) matrix and w is an (n,n) matrix defined as 


[P] = [p(*i) P(*2> ••• PUJ ] ' , 


w x (x) 


w 2 (x) 


w n (x) 


where x = x- Xj , and 


T 2 

p (*)= 1, X, X , 


p^(.*:)=— ^— = 0 , 1 , 2x, ... (m- \)x r 
dx 


with (m- 1) as the order of the basis function p(x) used in the MLS approximation. The 
weight functions w f (x) chosen are 


„ I" 1 -d?/R? 1 ' if d: <R: 

v< W = ) L ' ' J 


if d, > R, 


\ { d Uf (dA 4 

1-6 — +8 — -3 — if 0 <d:<R: 

Wi {X ) = \ R R R 


if df > Rf 


where dj = II x-x, II. The test function v,(x) in the MLPG weak form is chosen as 


v,: U) = 


1 ~d-/R' c 


if dj < R a 


if dj > R a 


Note that the lengths R, and R 0 in Eqs. (9) and (10) are user defined in the MLPG method. 
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Ill the current implementation a beam of length 4 1 is considered, as the choice of 
unit beam length l would mask numerical errors. Six models with 5, 9, 17, 33, 65, and 
129 nodes uniformly distributed along the length of the beam are considered. The model 
with 17 nodes is presented in Figure 1. The distance between the nodes (A//) in these 
models are 1.0, 0.5, 0.25, 0.125, 0.0625, and 0.03125 for the 5-, 9-, 17-, 33-, 65-, and 
129-node models, respectively. The ( R a / /) in the test functions (Eq. 10) in each of these 
six models is different and is chosen equal to (2A). The (R, / /) in Eq. (9) is chosen to be 
(Ri / / = 3.5) for the 5-, 9-, and 17-node models and (R, 11= 16A) for the 33-, 65-, and 
129-node models. Two types of basis functions, quadratic basis (7 , x, X 2 ) and cubic basis 
(/,.v,.v 2 ,,v“), are used. System matrices in Eq. (1) are developed with these parameters. 
The resulting system equations must be able to reproduce the constant, linear, and 
quadratic terms exactly when the quadratic basis is used, and additionally, the cubic term 
when the cubic basis is used. To evaluate the system matrices developed for the six 
models, two rigid body conditions and a constant-curvature condition were considered. 
These can be written as 

u(x) = Cq, 0 = — = 0 ; Rigid body translation, 

dx 

u(x) = CiX, 6 = Ci ; Rigid body rotation, and ( j ] ) 

2 

u(x)=c 2 X , 0 = 2c 2 x ; Constant - curvature, 

where Co, Cu and C 2 are arbitrary constants. The third condition in Eq. (11) corresponds to 

d 2 u 

the problem of a cantilever beam with a moment, M = El — = 2c 2 , applied at x = 4/. 

dx~ 

The problems described by Eq. (11) are simple test problems and should be reproduced 
exactly by the MLPG when quadratic or higher bases are used. 
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The { d } vectors that correspond to each of the conditions in Eq. (11) (and in the 


absence of any other loading) when used in Eq. (1) should result in a null right-hand 
vector if the K, (node ' ) is evaluated exactly, hi general, the product results in a residual 
{r {vector as 

K (n° de ) {d } = {r } (12) 

Each of the components of the vector { r } is nearly equal to machine zero if K/" od< J is 
evaluated accurately. To quantify the residual, an error norm of { r } is computed as 
ii hi rr^T 7 d3) 

111 \N dk= i 


where r k is the k' h component of the vector { r ) in Eq. (12) and Nd is the degrees of 
freedom in the model. 


Table 1; Error norm ||E'||i of the residuals for six models and for two basis functions 


Number of 
nodes in 
the model 

U~Ci 

U=C2X 

A, 

II 

Quadratic 

Basis 

Cubic 

Basis 

Quadratic 

Basis 

Cubic 

Basis 

Quadratic 

Basis 

Cubic 

Basis 


0.5040e-14 

0.1278e-12 

0.2099e-14 

0.4547e-13 

0.5733e-14 

0.9196e-13 


0.7515e-13 

0.1496e-ll 

0.2362e-13 

0.5514e-12 

0.3321e-13 

0.9747e-12 

17* 

0.2774e-10 

0.821 le-10 

0.1109e-10 

0.3067e-10 

0.1582e-10 

0.5352e-10 

33 

0.3609e-9 

0.1062e-5 

0.1266e-9 

0.4479e-6 

0.2587e-10 

0.9057e-6 

65 

0.1691e-6 

0.1435e-2 

0.7735e-7 

0.5855e-3 

0.1726e-6 

0.1193e-2 

129 

0.1796e-4 

0.5599e+0 

0.8154e-5 

0.2269e+0 

0.1794e-4 

0.4154e+0 


*R,/l = 3.5 


Table 1 presents the IIEIIi for the three conditions in Eq. (1 1) when the weight 
function in Eq. (9b) is used. (Similar results are obtained when weight function (9a) is 
used and hence these results are not presented here.) As seen from the table, the ILEIIi 
deteriorates with model refinement and for higher order basis. Closer examination of the 
residuals for each of the six models showed that the residuals were of machine accuracy 
for nodes near the origin while the residuals were largest at nodes farthest from the 
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origin. This observation was confirmed by running different cases with the origin at 
different locations along the length of the beam. Also, the residuals were largest for the 
models with the largest number of nodes. 

Closer scrutiny of computations showed that the numerical values of the shape 
functions for nodes that are systematically located about the center of the beam (for 
example, nodes 3 and 15, 2 and 16, and 1 and 17 in the 17-node model of Figure 1) are 
not exactly identical as expected. These differences increased with model refinement and 
when a higher basis was used. The error norm in Table 1 can be improved by using 
higher precision computations or inversion routines. However, a much simpler 
alternative to improve the accuracy is presented below. 

Local coordinate approach 

In the MLS interpolation, the basis functions are in terms of the global coordinate 
x. The [A] matrix thus formed using this basis is generally of the form (see Eq. 16, ref. 4) 

[A]=x{ w ^p-p‘ + m t(^p*-pI} (14) 

k = 1 


where x = x-Xj and M are the number of nodes in the domain of definition of node j for 
which the [A] matrix is being computed. (For convenience in presentation, the [A] 
matrices thus formed will be referred to as the global method.) As the order of the 
polynomial basis increases the conditioning of the [A] matrix deteriorates. For example, 
the matrix [A] will have terms like 7, A 2 , x 4 , x 6 on the diagonal for a cubic basis function. 
The [A] matrices for nodes near the origin and the [A] matrices for nodes farthest from 
the origin will be different. The conditioning is worse for [A] matrices for nodes farthest 
from the origin. This explains the differences in the error norms observed in Table 1. 
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The situation can be easily rectified if the MLS approximation is defined not in terms of a 
global basis, but rather in terms of a local basis. Figure 2 shows two identical shape 
functions, one centered at node /, and the other centered at node e. The global 
approximation for 

u(x) = p T (x)a(x) 

2 m — 1 

= rq + o 2 x + Q 3 x~ + + G m x 

can be rewritten in the neighborhood of node j, recognizing that x = Xj + E, where ^is a 
local coordinate measured from node j, as 

u(x) = a x + a 2 [xj + %)+ a 3 {x, +%)“+... 

= ( «t + a 2 xj + g 3 x 2 j j+ (a 2 + 2 a 3 Xj + . . .)% + {a 3 + . . .) f 

V , ^ (16) 
= b 1 +b 2 ^ + b 3 <^ +... 


where /;„ are the new undetermined coefficients in the MLS approximation. 

(A similar local coordinate transformation can be affected for node e in Figure 2 as 
x = x e + £.) The [A] matrix then is computed in a similar maimer as in Eq. (14) but with 


and 


p t (£> = 


pI(£> 


i, ^ ^ 2 , ... r _1 


0, 1, 14, 3<r, ... (m-l4 


(in — 2) 


(17) 


as A( ) = _!(). 

dx y ! ’ 


Local coordinate approach results 


Tire local coordmate approach is implemented in the evaluation of the shape 
functions and their derivatives for all the nodes in the six MLPG models of the beam. 
Table 2 compares the condition numbers of the [A] matrices at various locations on the 
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beam using global and local coordinate methods. The condition numbers are evaluated 
using routines available in NAPACK and the procedure outlined in references 5 and 6 . 
When the global coordinate method is used, the condition numbers of the [A] matrices 
for nodes farthest from the origin are much larger (suggesting poor conditioning) than the 
nodes closest to the origin. The conditioning numbers of the [A] matrices vastly improve 
when the local coordinate method is used, clearly demonstrating the advantages of the 
local coordinate method. 


Table 2: Comparison of the condition numbers of the [A] matrices at various 
locations on the beam using global and local coordinate methods 


Location 
on the 
beam (xJ4l) 

Number of nodes in the model 

5* 

9 * 

17* 

33 

65 

129 

Global Method Conditioning Number 


0.63 le+3 

0.106e+4 



0.267e+3 

0.189e+4 


0.23 le+5 

0.268e+5 

tWMSSM 


0.904e+6 

0.131e+8 

1.0 

0.914e+6 

0.77 le +6 

0.127e+7 

0.422e+8 

0.153e+10 

0.365e+ll 


Loca 

Method Conditioning Number 

0.0 

0.634e+3 

0.106e+4 

0.930e+3 

0.271e+3 

0.267e+3 

0.189e+4 

0.5 

0.478e+3 

0.496e+2 

0.41 le+2 

0.11 le +2 

0.153e+2 

0.141e+3 

1.0 

0.632e+3 

0.106e+4 

0.930e+3 

0.271e+3 

0.267e+3 

0.189e+4 


*Ri/l = 3.5 


Table 3: Error norm ||/f||i of the residuals computed with the local coordinate 
approach 


Number of 
nodes in 
the model 

U=Cj 

U=C2X 

A, 

II 

Quadratic 

Basis 

Cubic 

Basis 

Quadratic 

Basis 

Cubic 

Basis 

Quadratic 

Basis 

Cubic 

Basis 

5* 

0.1173e-14 

0.3500e-13 

0.2342e-15 

0.1201e-14 

0.3174e-14 

0.3853e-13 

9 * 

0.2521e-13 

0.4900e-13 

0.8357e-14 

0.1699e-13 

0.3659e-13 

0.4146e-13 

17* 

0.1392e-12 

0.2169e-12 

0.4764e-13 

0.1680e-12 

0.2126e-12 

0.8124e-12 

33 

0.4389e-12 

0.1390e-ll 

0.1876e-12 

0.5060e-12 

0.4084e-12 

0.2183e-ll 

65 

0.4196e-ll 

0.3890e-ll 

0.1142e-ll 

0.1879e-ll 

0.2548e-ll 

0.5930e-ll 

129 

0.4029e-10 

0.2778e-10 

0.1240e-10 

0.8191e-ll 

0.2400e-10 

0.2166e-10 


*Ri/l = 3.5 


The error nouns shown in Table 1 are recomputed and the results are presented in 
Table 3. As expected, all models and the quadratic and cubic basis functions produced 
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the error nouns close to machine accuracy, suggesting that the local coordinate approach 
produces accurate results compared to the global coordinate approach. 

Computational advantage of the local coordinate approach 

hi the conventional MLPG implementation, the [A] matrix is calculated and 
inverted at every node in the model. When using the local coordinate methodology with 
uniform nodal spacing, the shape functions are exactly identical for nodes whose Rj 
places the entire shape function in the interior of the domain of the problem. Hence, for 
those nodes the [A] matrices are identical. As such, considerable reduction in 
computational effort and cost can be achieved by the proposed local coordinate approach 
thus eliminating a perceived disadvantage of the MLPG method. 

Concluding remarks 

The MLPG method for beam problems (C 1 problems) showed that the solutions 
deteriorated as the number of nodes in the models were progressively increased. Closer 
examination revealed that the moving least squares (MLS) shape function calculations 
involved the computation of the [A] matrix and that this matrix became ill conditioned 
for nodes farthest from the origin. To overcome this difficulty a local coordinate 
approach for the MLS basis functions was proposed. The proposed approach restored the 
accuracy of the MLPG method for beam problems. 
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